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Abstract. 

Two models of anomalous difTusion of cosmic ray in the leaky-box approximation [T] are 
compared: one of them is based on the decoupled time-space Levy flights and the other on 
fractional walks with a finite free motion velocity. Distributions of first passage time and paths 
are computed and evolution of diffusion packets to equilibrium state is shown. Calculations 
demonstrate essential difference between the two models: the coupled scheme gives more realistic 
results. 



5^ 1. Introduction 

We continue investigation of problems arising in the fractional model of cosmic rays propagation 
in the Galaxy. Recall the situation. The first version of this model was proposed in f2j. It was 
based on the three-dimensional Levy-Feldheim flight process in an infinite homogeneous medium 
obeying the equation with a fractional power of Laplacian 



G(r,t) = 4(r)i(t), 



containing the "diffusion coefficient" D, which depends on the energy E and remains constant 
in the process of motion, and the fractional Laplacian given by its Fourier transform 

J e'^^-A)''/^f{r)dr = |k|" J e^^''f{r)dr = |k|°7(k), < a < 2. 

This equation and its solution expressed through the isotropic Levy-Feldheim distribu- 
tion *(")(x) 

G{r,t) = (L»t)-3/°^W(r(Dt)-3/°) 

were known to that time (see ^). The fractional power of the Laplacian has been explained by 
a fractal (friable) large scale structure of the Galaxy which causes an enhanced kind of diffusion 
(superdiffusion). The form of such a diffusion packet is described by the isotropic Levy-Feldheim 
distribution, and its width grows with time proportionally to t^/". When a = 2, the medium 
becomes homogeneous and anomalous diffusion reduces to the normal one with Gaussian profile 
and the width oc t^/^. As a result, we saw that this solution on the assumption of power 
type of energy dependence of diffusivity {D oc E^) and source spectrum {S oc E^'P) reveals 
an effect similar to the observed "knee" in primary spectra. As noted at the end of the cited 
work, the best fit of experimental data for H, He, CNO, Ne-Si and Fe-group were observed at 
a = 5/3 ~ 1.67, 6 = 0.25 and p = 2,9. We will refer to this model as the LU-model. 



2. Lagutin-Tyumentsev model 

Later, the model was modified by inserting a fractional time-derivative instead of the first-order 
one [3] and lowering the orders a of the fractional Laplacian to 0.3 [5]: 

G(r, t) = 5{Y)5p{t), 5p{t) = t~^/T{l - /?), a = 0.3, /3 = 0.8. 

For the sake of convenience, we will refer to this modification as the LT-model. 

The following case for such a choice was given. The value (3 = 0.8 was taken from the work [6] 
devoted to investigation how photospheric convective motions transport magnetic flux elements. 
Experimental data exhibit subdiffusion behaviour of solar magnetic elements. Observations 
of solar magnetic bright points analyzed in the cited work led to conclusion that the waiting 
time distribution density follows t~^~^, (3 = 0.61 it 0.09, during interval 0.3-22 minutes and 
then rapidly damps. Thus, we do not see here any reasons for application of the power law 
with f3 = 0.8 to description of propagation of cosmic rays through interstellar medium: these 
phenomena are quite different by space-time scales and even by nature. Moreover, the power- 
law behavior is observed at small times only and disappears in the long-time region. This 
corresponds to the value (3 = 1. 

The spatial exponent a was changed firstly from 1.67 to 1.00 in The authors wrote that 
"a comparison of simulated characteristics with experiment indicates that the fractal structure 
of ISM with the parameter a = 1 (Kraichnan spectrum of magnetic irregularities) gives local 
cosmic ray characteristics which are closest to the experiment". However, as far back as in 1999 
[8], see also was found, that the fractal dimension dp oi a medium does not coincide with 
the exponent a characterizing the free path distribution in this medium. This conclusion was 
supported by analytical calculations performed in [10]. They showed that a depends not only 
on dp but also on size of scattering objects (say, magnetic clouds). When it was recognized, 
Lagutin, Raikin and Tyumentsev [11] repeated the calculations and determined the spatial 
exponent a = 0.3 from the linearized relation 

a = 2 — dp 

with dp = 1.7 [11]. The latter number for the fractal dimension of the interstellar medium can 
be considered as a conventional value (see, for example, [12]). Link between a and dp is valid 
only in the limit of point fractal ([llj). Approximately, this formula could be used in the case 
if the interstellar magnetic field form small islands located at large distances of each other, but 
really it is not the case. Interstellar magnetic clouds affecting the cosmic ray transport have 
various sizes and may be close to each other. Thus, expecting shorter free paths, one should take 
essentially larger values of a. Looking at Figure 3 of the article [TT], one can see an ambiguous 
a{dp) dependence: the a is determined not only by the dp but also by ratio size/distance for 
inhomogeneities. When this ratio grows, the fractal becomes less transparent and the free path 
pdf falls more rapidly. We refer to the work [13], where the numerical simulation with the 
Erlykin-Wolfendale model gave a = 1.6 -i- 1.9 which, in the authors opinion, "is expected result, 
implying a fractal structure for interstellar medium". The choice is also supported by recent 
article [H]. Using the known cosmic ray spectrum and radial gradient in the vicinity of the solar 
system to define an energy density and comparing with the modeling results shoved that the 
best fit for the value of a is about 1.65 (possible fits range from 1.6-1.9, but not acceptable fit is 
found for a = 2, which would correspond to conventional diffusion). Our initial value a = 1.67 
belongs to this region. 

3. Bounded anomalous diffusion model 

In our works \15\ [TBI llZ] we show that random trajectories related to LT-model contains 
anomalously long rectilinear parts comparable with size of the Galaxy disc itself, and assumption 
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Figure 1. Examples of time-position trajectories for three models: normal diffusion (path length 
/ = 1 pc), anomalous diffusion (a = 1.67,(3 = 0.8, Dq = 2.4 • lO^^ pcl•Vyear°■^ D = DqR^-'^\ 

.0.3 , 



E = W GeV), and LT model (a = 0.3, /? = 0.8, Z^o = 4 • 10"^ pc"■Vyear"•^ D = DoR°-^\ 
E = 10^ GeV). The lower panels present schematically the diffusion packet spreading laws, 
shadowed regions correspond to regions with determining influence of speed finiteness. 



on instantaneous flights to such distances looks to say the least of it unphysical (see Figure [T]) . 
The exit from this situation lies in using the fractional material derivative operator as it described 
in the above-sited our works: this operator takes into account that cosmic ray particles propagate 
through interstellar medium with a finite speed, the fractional material operator used in our 
bounded anomalous diffusion model [TB] , where the following equation 



A^{oV^)G{r,t) = S,{r,t), < a < 1; 



ld_ 

V dt 



G{T,t) = S^{Y,t), l<a<2; 



for cosmic ray propagation was represented. Here, G{r,t) is the propagator, Sy{r,t) is a source 
function, the angle brackets denote averaging over directions fl of propagation, the operator 

-- + nvj G{r,t) 

is the fractional generalization of the material derivative. When a G (1,2), the equation reduces 
in the long time asymptotic region to the Levy-flight diffusion equation 

^ = -!?,(- A)"/2G(r,t), G(r,0)=<5(r), 



with diffusivity [19 



1 + w/v 

where w stands for he mean path covered by Levy-jumps per unit time and v is the free motion 
of particles. The solution of Eq. (2) for an unbounded fractal medium is expressed through the 
isotropic Levy-Feldheim distribution 

G{r,t) = (Z)„t)-3/°^H((i?^t)-i/"r), 1 < a < 2. 



4. Fractional Laplacian in a bounded domain 

In this work, we consider another aspect of cosmic ray propagation, provoked by the long- 
distant parts, namely the influence of boundaries on the fractional Laplacian. Indeed, the true 
Laplacian is a local operator, having the same form independently of presence or absence of 
boundaries. However, the fractional Laplacian is a non-local operator and for this reason it has 
a form depending on boundaries. In particular, the definition based on the Fourier transform 
can not be applied to the fractional Laplacian acting in a bounded medium. 

The statement of such a problem should be accompanied with a specification of the desired 
function values throughout an outer region. So, we have to return to the integral representation 
of the operator. The random flight interpretation can help in specifying the conditions but some 
subtle points such as distinction between first-passage and first-arrival times or between free and 
reflecting boundary conditions appear [2D]. In [20] have investigated the matrix representation 
of the one-dimensional fractional Laplacian and solved numerically in connection to the first- 
passage problem (the Levy-flights under absorbing boundary conditions) and to the long-ranged 
interfaces with no constraints at the ends (the free boundary conditions). 

Krepysheva et al ({21]) analyze the symmetric Levy flights restricted to a semi-infinite domain 
by a reflective barrier. They show that the introduction of the boundary condition induces a 
modification in the kernel of the nonlocal operator: 
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The operators — (— A)"/^ and —{—A)"^^ differ in the kernels, but the difference becomes small 
when a; + ^ is large. Nevertheless, omitting the term {x + we would get a decreasing 

integral with respect to x, whereas the total amount of the diffusing matter should be preserved. 

Rafeiro and Samko ([22]) have introduced a version of the fractional Laplacian for a bounded 
domain as a generalization of the Marchaud formula for one-dimensional fractional derivatives 
on an interval (a, 6), — 00 < a < 6 < 00, to the multidimensional case of functions defined on a 
region G CW^: 
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In other words, this is the Riesz fractional derivative of the zero continuation of /(x) from G to 
the whole space W^. 

Guan & Ma [23], investigating the reflected symmetric a-stable processes, gave the name 
regional fractional Laplacian to the limit 



-(-A)^/V(x) =lim C(a) 



G, |x-y|>e 



provided it exists. 

For more detail, the reader can be referred to the articles [Ml ESI [Ml 123 [Ml 129] • Better 
understanding of the Laplacian in a bounded domain can be achieved on the base of the non- 
local operator theory [30 l [3T 1 [32]. 

This short review is done in order to underline that in contradistinction to classical case, 
the fractional Laplacian A"/^ change its form in a bounded domain and cannot be determined 
by its Fourier transform — |k|° anymore. Consequently, all results obtained by Lagutin et all 
in 2001-2011 years relate to infinite unbounded fractal medium. One should say, that referring 
to the normal model with the use of Gaussian distribution in a bounded model can not justify 
the similar use of the stable distributions because their long tails may easily get the boundary 
surfaces which are inaccessible for the normal Gaussian process. 

5. Numerical simulation 
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Figure 2. The scheme of the model {h = 150 pc). 



Taking into account the above-mentioned difficulties with statement of boundary conditions 
in analytic or numerical approach, we perform direct Monte Carlo simulations to investigate the 
CR propagation in the framework of a fractal Galaxy model. The first problem we consider here 
is the escape time distribution for the Galactic disk. As can be concluded from [33] (see also 
[34] ) , the leading contribution in this process belongs to plane boundaries so the escape through 
the cylindrical part of the boundary can in the first approximation be neglected. Thus, we will 
simulate isotropic walk of particles in a fractal layer with two plane-parallel boundaries and the 
initial random point uniformly distributed between these boundaries (Figure [2]). Objects for 
study are escape time and escape path distributions in two models: LU (a = 1.67, /? = 1, 
V = c) and LT (a = 0.3, /3 = 0.8, v = oo). Results of Monte Carlo simulation are presented in 
Figure [31 The coefficient of anomalous diffusion for both models D = DqR^'^"^ , E = 10^ GeV, 
parameter Dq = 4 • 10~^ pc'^ '^/year'''^ in LT-model and Dq = 2.4 • 10~^ pc^/year*^'^ in LU- 
model. In the LT-model, distributions of free path lengths and waiting times have the form 
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Figure 3. Numerically calculated pdf of first passage time (left panel) and path (right panel) 
in frameworks of GS-1963, LU-2000 and LT-2004 models. 



of asymptotical power laws 

P{e > n ~ -^^f^, a>0, r^oo; P{r > t} ^ ^ ^ /? > t ^ oo. 
r(l — a) Til — p) 

We take Cq = 270 pc^^ and = 10~^ yr~^ (for E = 10^ GeV) and simulate ^ and r as random 
variables with pdf in the form of fractional exponents: 

|lnC/|^/" |lnC/|i/" 



Here S'(a) and S{I3) are one-sided stable variables simulated according to Kanter's algorithm 

d sin(a^C/2)[sin((l - a)^[/2)]i/"-^ 



S(a) 



[sin(7r;72)]V°[ln;73]i/ 



where C/i, C/2 and C/3 are variables uniformly distributed in (0, 1]. 

In the LU-model, waiting times can be simulated according to the exponential distribution 

P{r > t} = exp(-/it). 

We take /x = 10~^ yr"^- Path length are simulated according to the Pareto distribution 
P{? > r} = b r~" with 6 = 8- 10~^ pc~". These parameters are valid for E = 10^ GeV. 

Black dotted line corresponds to the ballistic motion from a source. Distance between planes 
is equal to 2h = 300 pc. Right panel of Figure 3 shows the numerically calculated pdf of first 
passage time. Blue solid lines are for LU-model and red dashed line is for LU model for the 
same parameters as in the left panel. 

The mean escape time can be calculating according to the following formula 

Tesc = £TSB + e(l - £)[rSB + Tbb] + e(l - efirSB + '^Tbb] + ■ ■ ■ = TSB + -TBB, 
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Figure 4. Escape time versus transparency. The instantaneous point source is situated on the 
middle plane. 
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Figure 5. Distribution density of the transverse coordinate z of the cosmic rays in the region 
with specularly reflecting boundaries z = ±h. The instantaneous point source is on the middle 
plane. 



where e is the transparency of boundaries, tsb is the mean passage time from source to boundary, 
and Tsb is the mean passage time from boundary to boundary. In Figure HI the dependences 
of escape time on transparency in the models under consideration are shown. Corresponding 
values of tsb and tbb ai^e indicated in the figure. In the LT-2004 model the mean escape time 
is infinite due to trapping times with asymptotically power law distributions. 

In our work ^17j, it has been shown by Monte Carlo simulation that the LT-model provides 
large anisotropy for cosmic rays propagated in infinite space from a single source. From Figure [3l 
we can see that even specularly reflecting boundaries can not change this situation. First 
passage time are distributed in very wide interval. Figure [5] confirms this reasoning. It shows 
distribution density of the transverse coordinate z for the case of random walk of a particle 
between two specularly reflecting boundaries with coordinates z = ±h. Random walk starts 
from the middle plane z = 0. The densities are calculated for several times. One can see that 



even for t = 10 yr stationary distribution is not established. For the LU-model, the uniform 
distribution of transverse coordinate takes place at time t ~ 5 • 10^ yr for parameters indicated 
above. 
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